Scenario definitions

  • Tree canopy

    • Population-based Scenario: AI: Increase by 10% in all zip codes

    • Targeted

      • Scenario AII1: Increase by 10% in zip codes in the lowest 1/5th of current TC cover (i.e. <=20th pctile)

      • Scenario AII2: Increase by 10% in zip codes in the highest 1/5th of the Social Vulnerability Index (i.e. >80th pctile)

      • Scenario AII3: Increase by 10% in zip codes in the highest 1/5th of hospitalization burden (i.e. >80th pctile)

    • Proportionate-universalism

      • Scenario AIII1: Increase by 10% for bottom 1/5th of current TC cover… down to 2% for top 1/5th

      • Scenario AIII2: Increase by 10% for top 1/5th of SVI … down to 2% for bottom 1/5th

      • Scenario AIII3: Increase by 10% for top 1/5th of hospitalization burden … down to 2% for bottom 1/5th

  • Impervious surface cover

    • Population-based: Scenario BI: Decrease by 10% in all zip codes

    • Targeted

      • Scenario BII1: Decrease by 10% in zip codes in the highest 1/5th of current imperv cover (i.e. >80th pctile)

      • Scenario BII2: Decrease by 10% in zip codes in the highest 1/5th of the Social Vulnerability Index (i.e. >80th pctile)

      • Scenario BII3: Decrease by 10% in zip codes in the highest 1/5th of hospitalization burden (i.e. >80th pctile)

    • Proportionate-universalism

      • Scenario BIII1: Decrease by 10% for top 1/5th of current imperv cover … down to 2% for bottom 1/5th

      • Scenario BIII2: Decrease by 10% for top 1/5th of SVI … down to 2% for bottom 1/5th

      • Scenario BIII3: Decrease by 10% for top 1/5th of hospitalization burden … down to 2% for bottom 1/5th

Data check

Number of zip codes with scenarios by scenario type

hosp_all_long %>% 
  filter(measure=="irr") %>% 
  group_by(scenario_intervention,scenario_type_7_abbrev) %>% 
  summarise(n_scenario=n()) %>% 
  ungroup()
## # A tibble: 14 × 3
##    scenario_intervention scenario_type_7_abbrev n_scenario
##    <chr>                 <chr>                       <int>
##  1 Imp                   PB                           1260
##  2 Imp                   PU 1                         1260
##  3 Imp                   PU 2                         1260
##  4 Imp                   PU 3                         1260
##  5 Imp                   T 1                          1260
##  6 Imp                   T 2                          1260
##  7 Imp                   T 3                          1260
##  8 Tree                  PB                           1260
##  9 Tree                  PU 1                         1260
## 10 Tree                  PU 2                         1260
## 11 Tree                  PU 3                         1260
## 12 Tree                  T 1                          1260
## 13 Tree                  T 2                          1260
## 14 Tree                  T 3                          1260

Implausible or unusual data Are there zip codes with ratios above 2? This implies the counterfactual scenario has twice as many cases as baseline

library(knitr)
n_by_scenario=hosp_all_long %>% 
    filter(measure=="irr") %>% 
    group_by(scenario_intervention,scenario_type_7_abbrev) %>% 
    summarise(n_scenario=n())

hosp_all_long %>% 
  filter(measure=="irr") %>%
  filter(value>2) %>% 
  group_by(scenario_intervention,
           scenario_type_3,scenario_sub_type, scenario_type_7_abbrev) %>% 
  summarise(n=n()) %>% 
  ungroup() %>% 
  left_join(n_by_scenario,by=c("scenario_intervention",
    "scenario_type_7_abbrev")) %>% 
  mutate(prop=n/n_scenario) %>% 
  kable(digits=4)
scenario_intervention scenario_type_3 scenario_sub_type scenario_type_7_abbrev n n_scenario prop
Imp Population-based NA PB 25 1260 0.0198
Imp Prop. Univ. Exposure PU 1 13 1260 0.0103
Imp Prop. Univ. Hosp. burden PU 3 10 1260 0.0079
Imp Prop. Univ. SVI PU 2 9 1260 0.0071
Imp Targeted Exposure T 1 7 1260 0.0056
Imp Targeted SVI T 2 3 1260 0.0024
Tree Population-based NA PB 9 1260 0.0071
Tree Prop. Univ. Exposure PU 1 2 1260 0.0016
Tree Prop. Univ. Hosp. burden PU 3 4 1260 0.0032
Tree Prop. Univ. SVI PU 2 3 1260 0.0024

Ratios below zero? This implies the counterfactual scenario yields a negative number (negative predicted hospitalizations)

hosp_all_long %>% 
  filter(measure=="irr") %>%
  filter(value<0) %>% 
  group_by(scenario_intervention,
           scenario_type_3,scenario_sub_type, scenario_type_7_abbrev) %>% 
  summarise(n=n()) %>% 
  ungroup() %>% 
    left_join(n_by_scenario,by=c("scenario_intervention",
    "scenario_type_7_abbrev")) %>% 
  mutate(prop=n/n_scenario)%>% 
  kable(digits=4)
scenario_intervention scenario_type_3 scenario_sub_type scenario_type_7_abbrev n n_scenario prop
Imp Population-based NA PB 20 1260 0.0159
Imp Prop. Univ. Exposure PU 1 15 1260 0.0119
Imp Prop. Univ. Hosp. burden PU 3 12 1260 0.0095
Imp Prop. Univ. SVI PU 2 9 1260 0.0071
Imp Targeted Exposure T 1 5 1260 0.0040
Imp Targeted SVI T 2 1 1260 0.0008
Tree Population-based NA PB 5 1260 0.0040
Tree Prop. Univ. Exposure PU 1 3 1260 0.0024
Tree Prop. Univ. Hosp. burden PU 3 3 1260 0.0024
Tree Prop. Univ. SVI PU 2 3 1260 0.0024

Plots

Histograms of IRR, IRD, and diff in n cases

Facet by type of intervention (impervious surfaces vs tree canopy) and by scenario type - Population-based, Proportionate Universalism, Targeted

facet_histogram_fun=function(df){
  df %>% 
      ggplot(aes(value))+
      geom_histogram()+
      facet_grid(
        rows=vars(scenario_type_7_abbrev),
        cols=vars(scenario_intervention)
      )
}

#IRR
hosp_all_long %>% 
  filter(measure=="irr") %>%
  #Remove ratios above 2
  filter(value<2) %>%
  #remove negative ratios
  filter(value>0) %>% 
  facet_histogram_fun()+
  xlab("IRR")

#IRD
hosp_all_long %>% 
  filter(measure=="ird") %>%
  filter(value<1e-6) %>% #Exclude outliers
  facet_histogram_fun()+
  xlab("IRD")

#n cases
hosp_all_long %>% 
  filter(measure=="n_cases_diff") %>%
  facet_histogram_fun()+
  xlab("n_cases_diff")

Bar charts

Note I took the median of each measure over all zip codes.

Median IRR

hosp_all_long %>% 
  filter(measure=="irr") %>% 
  group_by(measure, scenario) %>% 
  summarise(irr_med=median(value,na.rm=TRUE)) %>% 
  ungroup() %>% 
  left_join(lookup_scenario, by ="scenario") %>% 
  ggplot(aes(y=irr_med,x=scenario_type_7_abbrev))+
  geom_col()+
  labs(x="Scenario", y="Median IRR over ZCTA")+
  facet_grid(
    #Move facet down to x-axis
    #https://stackoverflow.com/questions/67519146/bar-plot-with-named-groups-on-x-axis-in-ggplot2
    cols=vars(scenario_intervention),
    scales="free_x",
    space="free_x",
    switch="x"
  )+
  theme(panel.spacing = unit(0, units = "cm"), # removes space between panels
        strip.placement = "outside", # moves the states down
        strip.background = element_rect(fill = "white") 
  )
## `summarise()` has grouped output by 'measure'. You can override using the
## `.groups` argument.

Median IRD

hosp_all_long %>% 
  filter(measure=="ird") %>% 
  group_by(measure, scenario) %>% 
  summarise(ird_med=median(value,na.rm=TRUE)) %>% 
  ungroup() %>% 
  left_join(lookup_scenario, by ="scenario") %>% 
  ggplot(aes(y=ird_med,x=scenario_type_7_abbrev))+
  geom_col()+
  labs(x="Scenario", y="Median IRD over ZCTA")+
  facet_grid(
    #Move facet down to x-axis
    #https://stackoverflow.com/questions/67519146/bar-plot-with-named-groups-on-x-axis-in-ggplot2
    cols=vars(scenario_intervention),
    scales="free_x",
    space="free_x",
    switch="x"
  )+
  theme(panel.spacing = unit(0, units = "cm"), # removes space between panels
        strip.placement = "outside", # moves the states down
        strip.background = element_rect(fill = "white") 
  )
## `summarise()` has grouped output by 'measure'. You can override using the
## `.groups` argument.

Median n_cases_diff

hosp_all_long %>% 
  filter(measure=="n_cases_diff") %>% 
  group_by(measure, scenario) %>% 
  summarise(n_cases_diff_med=median(value,na.rm=TRUE)) %>% 
  ungroup() %>% 
  left_join(lookup_scenario, by ="scenario") %>% 
  ggplot(aes(y=n_cases_diff_med,x=scenario_type_7_abbrev))+
  geom_col()+
  labs(x="Scenario", y="Median diff. in N cases over ZCTA")+
  facet_grid(
    #Move facet down to x-axis
    #https://stackoverflow.com/questions/67519146/bar-plot-with-named-groups-on-x-axis-in-ggplot2
    cols=vars(scenario_intervention),
    scales="free_x",
    space="free_x",
    switch="x"
  )+
  theme(panel.spacing = unit(0, units = "cm"), # removes space between panels
        strip.placement = "outside", # moves the states down
        strip.background = element_rect(fill = "white") 
  )
## `summarise()` has grouped output by 'measure'. You can override using the
## `.groups` argument.

Total n_cases_diff

hosp_all_long %>% 
  filter(measure=="n_cases_diff") %>% 
  group_by(measure, scenario) %>% 
  summarise(n_cases_diff_sum=sum(value,na.rm=TRUE)) %>% 
  ungroup() %>% 
  left_join(lookup_scenario, by ="scenario") %>% 
  ggplot(aes(y=n_cases_diff_sum,x=scenario_type_7_abbrev))+
  geom_col()+
  labs(x="Scenario", y="Sum diff. in N cases over ZCTA")+
  facet_grid(
    #Move facet down to x-axis
    #https://stackoverflow.com/questions/67519146/bar-plot-with-named-groups-on-x-axis-in-ggplot2
    cols=vars(scenario_intervention),
    scales="free_x",
    space="free_x",
    switch="x"
  )+
  theme(panel.spacing = unit(0, units = "cm"), # removes space between panels
        strip.placement = "outside", # moves the states down
        strip.background = element_rect(fill = "white") 
  )
## `summarise()` has grouped output by 'measure'. You can override using the
## `.groups` argument.

Interactive maps

IRR

Impervious surfaces

Population-based

Prop. Univ.